#ifndef THERMAL_FI_OPERATOR_CONNECTIONS_MPFA_2  
#define THERMAL_FI_OPERATOR_CONNECTIONS_MPFA_2

#include "rs_smesh_base.h"

#include "boundary_building.h"
#include "residual1D.h"
#include "residual2D.h"

namespace blue_sky{

	namespace mpfa{

template <typename strategy_t, bool is_w, bool is_g, bool is_o>
struct thermal_mpfa_base_impl : mpfa_base_impl<strategy_t, is_w, is_g, is_o>
{
	typedef mpfa_base_impl<strategy_t, is_w, is_g, is_o> base_t;		
	typedef typename strategy_t::csr_matrix_t           csr_matrix_t;
	typedef typename csr_matrix_t::i_vector_type		i_vector_type;
	typedef typename csr_matrix_t::fp_vector_type		fp_vector_type;
	typedef typename strategy_t::item_t					item_t;
	typedef typename i_vector_type::const_iterator		iterator;

	//constructor
	thermal_mpfa_base_impl (
	const thermal_fi_operator_impl <strategy_t, is_w, is_g, is_o> &fi_operator_,
		rhs_item_array_t &rhs_, rhs_item_array_t &reg_values_, 
		rhs_item_array_t &ss_values_)
	: base_t(fi_operator_, rhs_, reg_values_), temperature_(fi_operator_.temperature_),
	s_values_(ss_values_), ss_diag_(fi_operator_.ss_diag_), 
	ss_rows_(fi_operator_.ss_rows_ptr_), 
	cols_ind_(fi_operator_.ss_cols_ptr_),
	boundary_builder(fi_operator_.boundary_builder_)
	{
		NX = boundary_builder->dimens[0];
		NY = boundary_builder->dimens[1];

		shift = n_sec_vars - 1;
	}

	//return true if cell is on the boundary false otherwise
	bool is_on_boundary(index_t cell, coords boundary){
		i_vector_type boundaryList = *boundary_builder->get_result_map()[boundary];
		iterator begin = boundaryList.begin();
		iterator end	 = boundaryList.end();
		iterator result = std::find(begin, end, cell);
		return result != end;
	}

	/**
	* \brief  Fills s_rhs and diagonal part of ss_matrix                   
	* */
	void
	modify_boundary_sec_rhs (const double &dt)
	{						
		iterator it;									
		
		ResidualDiff1D<Residual1D<strategy_t>, strategy_t> residual1D(boundary_builder);
		ResidualDiff2D<Residual2D<strategy_t>, strategy_t> residual2D(boundary_builder);

		item_t value = 0;
		index_t cell;

		//left boundary
		i_vector_type boundaryX1 = *boundary_builder->get_result_map()[X1];
		i_vector_type boundaryX2 = *boundary_builder->get_result_map()[X2];
		i_vector_type boundaryY1 = *boundary_builder->get_result_map()[Y1];
		i_vector_type boundaryY2 = *boundary_builder->get_result_map()[Y2];

		for(it = boundaryX1.begin(); it != boundaryX1.end(); ++it){
			cell = *it;								
			if(NY > 1){									
				if( !is_on_boundary(cell, Y1) && !is_on_boundary(cell, Y2)){
					value = residual2D(cell, dt, temperature_[cell], 
							temperature_[cell + 1], temperature_[cell], 
							temperature_[cell + NX], temperature_[cell - NX], 
							dRdUij, dRdUipj, dRdUimj, dRdUijp, dRdUijm);										
				} else{
					value = 0;
					dRdUij = 0;
				}
				
			} else{
				value = residual1D(cell, dt, temperature_[cell], temperature_[cell + 1], 
					temperature_[cell], dRdUij, dRdUipj, dRdUimj);
			}								
			s_rhs_ [n_sec_vars * cell + shift] -=  value; 
			ss_diag_ [n_sec_vars * cell + shift] +=  dRdUij;
		}

		//right boundary
		for(it = boundaryX2.begin(); it != boundaryX2.end(); ++it){
			cell = *it;
			if(NY > 1){
				if(!is_on_boundary(cell, Y1) && !is_on_boundary(cell, Y2)){
					value = residual2D(cell, dt, temperature_[cell], 
					temperature_[cell], temperature_[cell - 1], 
					temperature_[cell + NX], temperature_[cell - NX], 
					dRdUij, dRdUipj, dRdUimj, dRdUijp, dRdUijm);
				} else{
					value = 0;
					dRdUij = 0;
				}
			} else{
				value = residual1D(cell, dt, temperature_[cell], temperature_[cell], 
					temperature_[cell - 1], 
					dRdUij, dRdUipj, dRdUimj);									 
			}
			s_rhs_ [n_sec_vars * cell + shift] -=  value;
			ss_diag_ [n_sec_vars * cell + shift] +=  dRdUij;
		} 	

		//inner cells
		i_vector_type non_boundary_cells = boundary_builder->get_non_boundary_cells();
		for (it = non_boundary_cells.begin(); it < non_boundary_cells.end(); ++it) {
			cell = *it;
			if(NY > 1){
				value = residual2D(cell, dt, temperature_[cell], temperature_[cell + 1],
					temperature_[cell - 1], temperature_[cell + NX], temperature_[cell - NX], 
					dRdUij, dRdUipj, dRdUimj, dRdUijp, dRdUijm);
			} else {
				value = residual1D(cell, dt, temperature_[cell], 
					temperature_[cell + 1], temperature_[cell - 1], 
					dRdUij, dRdUipj, dRdUimj);									
			}
			s_rhs_[n_sec_vars * cell + shift] -= value;
			ss_diag_ [n_sec_vars * cell + shift] +=  dRdUij;
		}

		
		if(NY > 1){
			//upper boundary
			for(it = boundaryY1.begin(); it != boundaryY1.end(); ++it){
				cell = *it;
				item_t left = cell - 1;
				item_t right = cell + 1;
				if(is_on_boundary(cell,X1)){
					left = cell;
				} else if(is_on_boundary(cell,X2)){
					right = cell;
				}

				value = residual2D(cell, dt, temperature_[cell], temperature_[right], 
					temperature_[left], 
					temperature_[cell + NX], temperature_[cell], 
					dRdUij, dRdUipj, dRdUimj, dRdUijp, dRdUijm);
				s_rhs_ [n_sec_vars * cell + shift] -=  value; 
				ss_diag_ [n_sec_vars * cell + shift] +=  dRdUij;
			}

			//bottom boundary
			for(it = boundaryY2.begin(); it != boundaryY2.end(); ++it){
				cell = *it;
				item_t left = cell - 1;
				item_t right = cell + 1;
				if(is_on_boundary(cell,X1)){
					left = cell;
				} else if(is_on_boundary(cell,X2)){
					right = cell;
				}

				value = residual2D(cell, dt, temperature_[cell], temperature_[right],
					temperature_[left], 
					temperature_[cell], temperature_[cell - NX],
					dRdUij, dRdUipj, dRdUimj, dRdUijp, dRdUijm);
				s_rhs_ [n_sec_vars * cell + shift] -=  value; 
				ss_diag_ [n_sec_vars * cell + shift] +=  dRdUij;
			}
		}

	}
		
	//fill non-diagonal part of ss_matrix for cell
	void calculate_element(const index_t cell, const index_t cell_ind, 
		const index_t cell_col_ind, index_t other_col_ind, 
		const double &dt,  index_t other_cell, index_t NX, index_t NY){
		
			ResidualDiff1D<Residual1D<strategy_t>, strategy_t> residual1D(boundary_builder);
			ResidualDiff2D<Residual2D<strategy_t>, strategy_t> residual2D(boundary_builder);
			const item_t zero = 0.0;

			if (other_col_ind == cell_col_ind + 1 && ss_rows_[cell + 1] > cell_ind + other_cell){
				if(NY > 1){
					residual2D(cell, dt, temperature_[cell], temperature_[cell + 1], zero, 
						zero, zero, dRdUij, dRdUipj, dRdUimj, dRdUijp, dRdUijm);
				} else {										
					residual1D(cell, dt, temperature_[cell], temperature_[cell + 1], zero,  
						dRdUij, dRdUipj, dRdUimj);											
				}	
				ss_values_ [n_sec_vars * (cell_ind + other_cell) + shift] +=  dRdUipj;	
			} else if (other_col_ind == cell_col_ind - 1){
				if(NY > 1){	
					residual2D(cell, dt, temperature_[cell], zero, temperature_[cell - 1], 
						zero, zero, dRdUij, dRdUipj, dRdUimj, dRdUijp, dRdUijm);
				} else {
					residual1D(cell, dt, temperature_[cell], zero, temperature_[cell - 1], 
						dRdUij, dRdUipj, dRdUimj);									
				}
				ss_values_ [n_sec_vars * (cell_ind + other_cell) + shift] +=  dRdUimj;
			} else if (other_col_ind == cell_col_ind + NX && 
				ss_rows_[cell + 1] > cell_ind + other_cell){
				if(NY > 1){
					residual2D(cell, dt, temperature_[cell], zero, zero, 
						temperature_[cell + NX], zero, dRdUij, dRdUipj, dRdUimj, 
						dRdUijp, dRdUijm);
					ss_values_ [n_sec_vars * (cell_ind + other_cell) + shift] +=  dRdUijp;	
				}
			} else if (other_col_ind == cell_col_ind - NX){
				if(NY > 1){
					residual2D(cell, dt, temperature_[cell], zero, zero, 
						zero, temperature_[cell - NX], dRdUij, dRdUipj, dRdUimj, 
						dRdUijp, dRdUijm);
					ss_values_ [n_sec_vars * (cell_ind + other_cell) + shift] +=  dRdUijm;
				}
			}
	}
	
	//fill non-diagonal part of ss_matrix for all cells
	void modify_ss_matrix (const double &dt)
	{	
		index_t total = boundary_builder->number_of_cells;
		index_t size = cols_ind_.size();
		index_t value;

		for (index_t cell = 0; cell < total; ++cell){
			const item_t cell_ind = ss_rows_[cell];
			const index_t cell_col_ind = cols_ind_[cell_ind];
			if(cell_ind + 1 < size){
				value = cols_ind_[cell_ind + 1];
				calculate_element(cell, cell_ind, cell_col_ind, value, dt, 1, NX, NY);
			}
			if(cell_ind + 2 < size){
				value = cols_ind_[cell_ind + 2];
				calculate_element(cell, cell_ind, cell_col_ind, value, dt, 2, NX, NY);
			}
			if(NY > 1){
				if(cell + 3 < size){
					value = cols_ind_[cell_ind + 3];
					calculate_element(cell, cell_ind, cell_col_ind, value, dt, 3, NX, NY);						
				}
				if(cell + 4 < size){
					value = cols_ind_[cell_ind + 4];
					calculate_element(cell, cell_ind, cell_col_ind, value, dt, 4, NX, NY);
				}
			}
		}
	}

	index_t NX, NY;
	index_t shift;
	smart_ptr<build_boundary<strategy_t>, true> boundary_builder;
	item_t dRdUij, dRdUipj, dRdUimj, dRdUijp, dRdUijm;
	const item_array_t              &temperature_;
	rhs_item_array_t                &ss_values_;
	rhs_item_array_t				&ss_diag_;
	const index_array_t				&ss_rows_;
	const index_array_t				&cols_ind_;
};

		template <typename fi_operator_t, typename mpfa_t>
		void
		thermal_block_connections_mpfa (typename mpfa_t::item_t dt, 
		fi_operator_t &fi_operator, mpfa_t mpfa_)
		{
				block_connections_mpfa(dt, fi_operator, mpfa_);                        

				typedef typename mpfa_t::index_t      index_t;
				typedef typename mpfa_t::rhs_item_t   rhs_item_t;
				typedef typename mpfa_t::item_t		  item_t;
				typedef boost::array <index_t, 3>     index_point3d_t;
				typedef mesh_rs<fi_operator_t::strategy_t> mesh_rs_t;
				typedef smart_ptr<mesh_rs_t, false>   sp_mesh_rs_t;

				const rhs_item_t *geom = &fi_operator.jmatrix_->geom_matrix->get_values()[0];
				
				mpfa_.m_mem = &fi_operator.m_array_[0];
				mpfa_.p_mem = &fi_operator.p_array_[0];
				
				data_dimens dims =  fi_operator.dimens_;
				index_t nx = dims.nx;
				index_t ny = dims.ny;
				mpfa_.modify_boundary_sec_rhs(dt);
				mpfa_.modify_ss_matrix(dt);
		}
			
	} // namespace mpfa


/**
* \brief  Calculates flux part of Jacobian
* \param  dt
* \return True
* */
template <typename strategy_t, bool is_w, bool is_g, bool is_o>
bool
thermal_fi_operator_impl <strategy_t, is_w, is_g, is_o>::block_connections_mpfa (
	const item_t &dt)
{
		mpfa::thermal_block_connections_mpfa (dt, *this,
			mpfa::thermal_mpfa_base_impl <strategy_t, is_w, is_g, is_o> (
			*this, flux_rhs_, reg_values_, ss_values_)
		);
		return false;
}
}

#endif //THERMAL_FI_OPERATOR_CONNECTIONS_MPFA_2